function [beta,lambda,se_beta,t_beta,lambda_se,a,r2,r2_2step]...
    =cross_section(test,factor)
% This function: two-step cross-sectional asset pricing test
% Output
% beta N2*N1 the beta estimate for each asset's first-stage regression
% lambda: 
% se_beta N2*N1 the standard error of beta estimate for each asset's first
% stage regression
% a N1*1 the column of constant 
% t_beta N2*N1 the tstats of beta estimates for each asset's first stage
% regression
% r2 N1*1 the R2 of each asset's first stage regression
% lambda N2*1 the price of risk for each factor
% lambda_se N2*1 the price of risk for each factor
% r2_2step the 2nd stage R2

[T,N1] = size(test); % number of test assets
[~,N2] = size(factor); % number of factors

% 1. sample selection: get rid of observations with all nan's in the
% cross-section
nan_test = sum(isnan(test),2); % each period how many nan
mm = max(find(nan_test == N1)); % the most recent period when obs of nan equals the size of the cross section
if isempty(mm) == 0 % if mm is not empty (exist observations with the whole cross section being nan)
factor = factor(mm+1:end,:);
test = test(mm+1:end,:);
T = T-mm;
end
% Note: the step above requires continuous observations --- part of the
% sample might be missing in early days, but it cannot allow for
% observations with gaps

% First step, make rooms
a = zeros(3,N1);
beta = zeros(N2,N1);
se_beta = zeros(N2,N1);
t_beta = zeros(N2,N1);
r2 = zeros(N1,1);
eps = zeros(T,N1);

for i=1:N1
    result=fitlm(factor,test(:,i));
    a(1,i)=result.Coefficients.Estimate(1,1);
    a(2,i)=result.Coefficients.SE(1,1);
    a(3,i)=result.Coefficients.tStat(1,1);
    beta(:,i)=result.Coefficients.Estimate(2:end,1);
    se_beta(:,i)=result.Coefficients.SE(2:end,1);
    t_beta(:,i)=result.Coefficients.tStat(2:end,1);
    eps(:,i)=result.Residuals.Raw;
    r2(i,1)=result.Rsquared.Ordinary;
end
exp_ret=nanmean(test,1)'; % average of asset returns using the unbalanced sample
% N1*1 matrix
beta = beta'; % N1*N2 matrix

% Without intercept in the second stage
result_2=fitlm(beta,exp_ret,'Intercept',false);
lambda=result_2.Coefficients.Estimate;
yfit = beta*lambda;
r2_2step = sum(yfit.*yfit)/sum(exp_ret.*exp_ret);

% Shanken's adjustment standard error
demean_f=factor-nanmean(factor,1);
sigmae = nw_balance(eps,0); % residual variance-covariance matrix, only using balanced sample
sigmaf = nw_balance(demean_f,0); % residual variance-covariance matrix, only using balanced sample

% Variance: Shanken's formula
var_lambda=((beta'*beta)\beta'*sigmae*beta/(beta'*beta)*(1+lambda'/sigmaf*lambda)+sigmaf)/T; % see Cochrane p240
lambda_se=sqrt(diag(var_lambda));
beta = beta';
end